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Abstract 



We provide an expression quantitatively describing the specific heat of the Ising model on the 
simple-cubic lattice in the critical region. This expression is based on finite-size scaling of numerical 
results obtained by means of a Monte Carlo method. It agrees satisfactorily with series expan- 
sions and with a set of experimental results. Our results include a determination of the universal 



amplitude ratio of the specific-heat divergences at both sides of the critical point. 
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I. INTRODUCTION 



Though real magnetic systems were supposed to be Heisenberg-like, the Ising model was 
originally introduced l| as a simplified model of magnetic ordering, because its relative 
simplicity offers better possibilities for a theoretical analysis. In later years, it was found, 
however, that Ising-like magnetic systems do exist. This is because real systems consist of 
spins embedded in a crystal lattice, and the resulting anisotropy field due to the neighboring 
charges may lift the 0(3) symmetry of an unperturbed spin. Depending on the character of 
the perturbation, the spin may have an 'easy axis' or an 'easy plane'. Here we consider the 
former case, which leads to Ising-like behavior. 

In many cases, the perturbation is relatively small and the system will approximately 
behave Heisenberg-like, except near an ordering transition where the paramagnetic state 
transforms into a lone-range ordered one. Near the transition, crossover occurs to 

Ising-like behavior. The critical singularities are then described by the Ising set of critical 
exponents. In some other cases, the perturbation due to the crystal field is so strong that 
the magnetic spins assume a true Ising character. This situation occurs when the ionic 
angular momentum S is described by a spin quantum number S > |, and the crystal field 
lifts the degeneracy of the S z eigenstates such that the S z = ±S doublet is lowest in energy, 
with the higher levels so far away that they play no role, even in the presence of exchange 
interactions between neighboring spins. Then the low-lying doublet can be described by an 
effective spin- 1/2 Ising Hamiltonian. This situation is known to occur for the Co 2+ ion in 
a tetrahedral coordination. It occurs also for some rare-earth ions like Dy 3+ and Yb 3+ in a 
sufficiently strong crystal field, with the provision that here the magnetic moments are due 
to spin as well as orbital angular momentum, and should thus be denoted J instead of S. 

If such ions are embedded in a crystal structure for which theoretical predictions for the 
thermodynamical properties such as the specific heat exist, comparison with experiments 
may be possible ji], f| . Such comparisons were made for dysprosium phosphate jf], Q] and for 
some alkali cobalt halides (si, [s| . These systems were found to behave, at least approximately, 
as the Ising models on the diamond lattice and the simple-cubic lattice respectively. 

The best way to obtain theoretical results for the thermodynamic properties of these 
models would obviously be an exact solution, but this is known to be a very difficult problem. 
It is thus noteworthy that it was 
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solution was found for the three-dimensional Ising model. However. Perk [11] and Wu ct 
al. [12I pointed out that Zhang's result for the free energy and the underlying arguments 
are flawed. Here it may be added that Zhang's result for the critical point of the simple- 



cubic Ising mode 
estimates 
accuracies 



13 



13 



s is not compatible with independent and mutually consistent numerical 
14]. The difference with Zhang's result exceeds the estimated numerical 
141 ] by several orders of magnitude. 
In the absence of an exact solution, one may still resort to approximations. At temper- 
atures sufficiently far above and below the critical point, excellent approximations exist in 
the form of series expansion of the partition function or the free energy, such as given in 



Refs. 
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and model on the simple-cubic lattice. In the critical region, the series of 

a finite length become inaccurate, and a method to extrapolate these series on the basis of a 
critical scaling assumption, such as used by Butera and Comi [17|, is needed. In the case of 
Rb 3 CoCl 5 (rubidium cobalt chloride) Q] the required theoretical prediction for the specific 
heat near criticality was also obtained this way. A similar analysis has been performed for 
the specific heat of DyP04 (dysprosium phosphate) p, LZI], which was instead compared with 
series expansions for the diamond lattice. However, these specific-heat analyses were con- 



known, for instance, 
was not included. 



ducted at a time that the value of the critical exponent a was not wel . 
a was set to zero in Ref. |9|. Moreover, Wegner's correction to scaling ] 

In order to obtain accurate predictions for the heat capacity in the critical region, one 
may apply Monte Carlo simulations. Cluster simulation methods 19|, [20], which strongly 
reduce critical slowing down, allow statistically accurate simulations in the critical region. 
Extrapolation of the finite-size simulation data to the thermodynamic limit is possible if the 
simulations cover a range of finite sizes exceeding the correlation length. Whereas this still 
excludes, as a result of the divergence of the correlation length, a narrow temperature range 
about the critical point, one may attempt to describe the extrapolated data by means of 
a scaling formula. The present work reports our efforts along this line for the case of the 
energy and the specific heat of the Ising model on the simple-cubic lattice. 

In Sec. Hi] we describe our Monte Carlo simulations, and the extrapolation to infinite 
system size. The derivation of scaling formulas for the energy and the specific heat, and the 
data analysis in terms of these formulas, are presented in Sec. IIHI Section [IV] discusses the 
numerical accuracies, provides comparisons with results from series expansions and with a 
set of experimental results, and ends with a few concluding remarks. 
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II. NUMERICAL TECHNIQUE 



The reduced Hamiltonian (Hamiltonian divided by kT) of the Ising model is denoted 

H{K) = -K^SiSi (1) 

M 

where the indices i and j label nearest-neighbor lattice sites on the simple-cubic lattice. The 
sum is on all nearest-neighbor pairs, and the spins Sk can assume values ±1. The coupling 
is defined by K = J/kT where J is minus the energy of a pair of parallel nearest-neighbor 
spins, k the Boltzmann constant, and T the temperature. The canonical reduced free energy 
density / is equal to 

/ = llnZ, Z = £V WW (2) 

{8} 

where Z is the partition function, iV the number of spins, and the sum is on all spin 
configurations {S}. The energy E and the specific heat C per particle, as expressed in 
dimensionless units, follow from the derivatives of / to K: 

E E df C = 2 fl 

J kTK dK" k dK 2 ' U 



A. Monte Carlo calculations 



Substitution of Eqs. (j2J) and (OQ) in Eqs. ([3]) leads to 

f - (4) 

and 

f = ^W 2 )-(H) 2 ), (5) 
where the ensemble averages (x), which are defined as 

{S} 

can be sampled directly using importance sampling. 

The simulations involved the sampling of the energy, as well as its square, for L x L x L 
Ising systems on simple-cubic lattices, with periodic boundary conditions. The system sizes 
were chosen as powers of 2 in the range 4 < L < 128, and in addition as L = 6 and 12. 
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About 10 7 samples were taken for L < 16, 2 x 10 6 for L = 32, 3 x 10 5 for L = 64, and 
5 x 10 4 for L = 128. Each sample was preceded by a number or Wolff cluster steps and/or 
Metropolis swee ps, depending on the value of K in comparison with the critical coupling 
K c m 0.2216546 [14j. For K « K c , Wolff clusters tend to be very small and only Metropolis 
sweeps were applied, and for K > K c only Wolff cluster steps. In the intermediate range, a 
few Metropolis sweeps were supplemented with a number of Wolff cluster steps. The number 
of Wolff clusters was chosen roughly equal to the inverse of the relative Wolff cluster size. 
The coupling K was given some 50 different values chosen to cover a wide range about the 
critical point. 



B. Extrapolation 

The analysis of the numerical finite-size data was done on the basis of well-documented 
finite-size scaling methods [21 ] . For non-critical systems with sizes L exceeding the correla- 
tion length, the data for the energy should approximately behave as 

E(K, L) = E(K, oo) + a(K)e- L/ ^ K) + ■■■ (7) 

from which the extrapolated energy E(K, oo) was obtained by means of a least-squares 
analysis. A small-system-size cutoff was applied when necessary to obtain a satisfactory 
residual x 2 . This cutoff varied between L = 6 far away from the critical point, and L = 32 
at a distance \K — K c \ « 0.005 from the critical point. No reliable extrapolations were 
obtained for \K — K c \ less than a few times 10~ 3 , with the exception of K = K c , where one 
expects that the finite-size-dependent energy converges as a power of L, which again enables 
extrapolation to L = oo. Typical estimated accuracies of the extrapolated results for E/J 
are in the order of 10~ 5 . 

The same extrapolation procedure was applied to the finite-size data for the specific heat 
with \K — K c \ > 0.005. Typical accuracies of the extrapolated results for C/k are estimated 
as at most a few times 10 -4 for K < 0.2 and K > 0.25, up to a few times 10 -3 in the vicinity 
K c . The extrapolated data are listed in the Appendix. 
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III. SCALING AND LEAST-SQUARES ANALYSIS 



A. Derivation from renormalization theory 

The analysis of the extrapolated data was done on the basis of scaling as derived from 
renormalization theory. The relevant equations follow from the assumptions that the picture 
described in the following paragraph is valid. 

The free energy density f(Ti,T 2 , ■ ■ ■) of the infinite system, expressed as a function of 
thermodynamic parameters Tj (J = 1, 2, • • •), can be written as the sum of an analytic 
part /.pi, T„ • • ■ ) and a singula, part /,. The singular part can be written f,( tl , t 2 , ■ ■ ■ ) as 
a function of Wegner's [22] nonlinear scaling fields tj, which are analytic functions of the Tj 
in a neighborhood of a critical point under investigation. Thus 

f(T u T 2 , ■ ■ ■ ) = f a (T u T 2 , ■ ■ • ) + f.ttuh, ■■■) (8) 

The singular part satisfies the scaling equation as implied by the renormalization theory. A 
rescaling of the linear dimensions by a factor b thus leads to 

f s {tiM,---) = h- d f s {WH u wn 2 ,---) (9) 

where d is the dimensionality and the yj are the renormalization exponents associated with 
the scaling fields tj, with the temperature exponent y\ positive, and the other exponents 
negative. The choice b = thus yields 

h, ■ ■ ■ ) = \t^f 8 {±l t Ihl-n^h, ■■■) (10) 

where ±1 has the sign of t±. Furthermore, f s (xi, x 2 , x$, ■ ■ ■ ) is an analytic function in a 
neighborhood of the X\ — 1, x 2 — 0, x% — 0, • • • . 

On the basis of this set of assumptions, we may Taylor expand the free energy in powers 
of the arguments Tj and tj, and then expand the tj's in the I}'s, resulting in an expres- 
sion depending only on the physical temperature fields, but with expansion coefficients that 
remain to be determined. We follow this procedure, restricting number of scaling fields in 
the expansion of Eq. (flOj) to two, namely the temperature field t = t\ and the irrelevant 
field u = t 2 . The corresponding exponents are denoted y t and y u respectively. The temper- 
ature exponent yt determines the leading singularity in the temperature-induced ordering 
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transition, while the irrelevant exponent y u generates Wegner's correction to scaling 18]. 
Expansion of the right-hand side of Eq. fflOl) thus yields 

f s (t,u) = |t|^^(j!)-V s 0j (±l,0)|t]- jW ^, (11) 

3 

where is the j'th derivative of f s to its second argument. The scaling fields are expanded 
as analytic power series in the temperature-like parameter to, defined by 

t = AK/K , AK = K-K C . (12) 

The analytic part of the free energy f a can be expanded directly in powers of AK. The 
resulting expansion of the total free energy density in powers of AK and t can be expressed 
in K, the only variable physical temperature parameter in our problem, as given by the 
Hamiltonian (Tj[|). Differentiation of the resulting expansion of the free energy density to K 
yields the dimensionless energy E/J. For d = 3 dimensions, the leading terms are included 
in 

-E(K, oo)/ J = 

V eJAKY + ffl a±\tf- yt)/yt + b ± u\tf- yt ~ y ^/ yt + p ± u 2 \t\ i3 ~ yt ~ 2y ^ /yt + ■■■ . (13) 
/ — ' dK 

3=0,1,- 

where we have included the first three terms in the sum on j in Eq. (JTT1) . and u is an analytic 
function of to related to u by 

^^(±1,0) §^u = b ± u (14) 
Vt dK 

The dimensionless specific heat C/k of the model ([1]) satisfies 

C(K, oo) _ ^ d 2 f(K, oo) _ K 2 dE 

k " dK* ~ TdK [Lb) 

and its expansion thus follows by differentiation of the energy, Eq. (fl3l) . This leads to 

3-y«-y« M^ ± | t |(3-2^ U )M + *i 6 itio-w-wuj/tt + . . . ( i 6) 
y t cm cm 

The parameters t and u, and their derivatives as they appear in Eqs. (fl3|) and (fl6|) . are 
expanded in powers of to as 



3 b ' 



= ±fl E iU - iHC 2 t § E ^ ■ (it) 

i=2,3,- i=i,2,- 
where ± stands for the sign of t, =F for its opposite, and 

«= X] «i<o» % = % J2 frA' 1 - (18) 

i=o,i.- '" i=i,2,- 

The scales of t and u are determined by setting it?! = it = 1. 



B. Fits 



Whereas Eq. (TTTT) includes, in principle, infinitely many terms, for numerical work it is 
necessary to truncate the expansion of / s , as well as those of f a and the scaling fields, at 
a finite order. Expression (1131) for the energy already contains the implicit simplification 
that there is only one irrelevant field, and that the expansion of f s (±l,x) is truncated at 
second order. Moreover, higher orders in the expansion of the temperature derivative of the 
irrelevant field were neglected. We shall reconsider these simplifications in Sec. IIV A[ No 
further simplifications were made in the derivation of Eq. ffT6]) from Eq. ffT3]) . 

Many attempts were made to fit Eqs. ( Tl3l) and (fl6l) to the numerical data, using different 
ranges of K, and different sets of parameters as determined by the orders at which the 
expansions were truncated. The unknown parameters in each set were determined by means 
of a Levenberg-Marquardt nonlinear least-squares analysis. Since Eqs. (IT31 and (1TB1) depend 
on the same parameters, the data for the energy and the heat capacity were simultaneously 
fitted by one set of parameters. 

A fit was considered satisfactory if it met three criteria: first, the residual x 2 has to be 
compatible with the number of degrees of freedom; second, there should be sufficiently large 
ranges of overlap with the accurate predictions from the low- and high-temperature series 
expansions; and third, at least the amplitudes of the leading terms in the fit formulas should 
be reasonably stable under variations of the i^-interval and of the number of correction 
terms in the temperature field and the analytic background. In Table [J we list the smallest 
satisfactory set of parameters thus obtained. We skipped the ellipses in Eqs. ( 1T31) and (TTB1) . 
and included terms up to order j = 4 in the expansion of t, up to j = 2 in that of u, and up to 
j = 5 in the analytic parts expressed by the first sums in Eqs. ffl3|) and ffl6|) . The residual of 
this fit was x 2 — 53.5, to be compared with the number of degrees of freedom df = 84. Since 
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possible correlations between specific heat and energy data could influence the estimation of 
the errors in the fitted parameters, we have analyzed the correlations between the deviations 
of the energy and of the specific heat with respect to the fit formula. We find a correlation 
coefficient of -0.066 which is not significant, and does not provide a reason to reconsider our 
error estimates. 

During the least-squares analysis, we found that some parameter values changed signifi- 
cantly when the A"-interval and/or the numbers of parameters in the expansions of t and of 
the analytic background were varied. Such shifts were sometimes comparable to the error 
margins as estimated from statistics based on the accuracy of the Monte Carlo results. This 
applied in particular to those of the Wj and the with j > 2. In this respect the amplitudes 
a + , a_, e , t\ and, to some extent, b + and b_ were better behaved. The error estimates 
listed in Tabled take into account the variation of the parameter values between these fits. 



IV. DISCUSSION 



A. Choice of parameters and their error margins 

Equation ffTTl) and the fits of E and C use only one irrelevant field, while, according to 
Newman and Riedel 23| . corrections to scaling could also arise from a second irrelevant 
field u' with exponent y u > ~ 2y u . We note that corrections generated in first order of u' 
would thus, in the present context, be practically indistinguishable from those generated 
in second order by u. For this reason, we have not included a separate term containing u'. 
Furthermore, the energy, Eq. f|T3|) . neglects a contribution due to the possible A"-dependence 
of the irrelevant field. Such a term behaves as \t\( 3 ~ yu ^ Vt and is thus a factor \t\ smaller than 
the leading correction. The third-order correction in u, which is also neglected, has nearly 
the same exponent. 

Another correction that was neglected is one with an integer exponent y" = —2, associated 
with the discreteness of the cubic lattice. The presence of such corrections could modify 
the higher-order correction amplitudes given in Table (TJ but the x 2 criterion did not yield 
indications that a term with y" = —2 should be included. 

Some insight in the relative importance of the corrections due to different orders of the 
irrelevant field can be obtained by comparing the fit including the second order of u, as given 
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TABLE I: Values of the parameters in the fit according to Eqs. (|13p and (|16p to Monte Carlo 
data in the interval 0.15 < K < 0.60. The error estimates given in the last column are not only 
based on statistics, but also on the variations of the parameter values due to changes of the fit 
interval and the number of parameters. In two cases the estimated error exceeds the parameter 
value and no error is quoted. While these values have no physical meaning, they are still useful for 
the evaluation of the specific heat and the energy. The values of yt, y u , and K c were taken from 
Ref. 14. 



parameter 


value 


error margin 


W\ 


1 


fixed 




0.662300 


0.06 


w 3 


0.160415 


0.09 




0.008397 




uo 


1 


fixed 


U\ 


-2.673700 


0.6 


CL- 


1.466642 


0.016 


a+ 


2.758572 


0.012 


b- 


0.923100 


0.2 


b+ 


-2.694381 


0.4 


P- 


-1.440041 


0.4 


P+ 


-2.345305 


0.8 


eo 


0.990604 


0.000004 


ei 


-27.847250 


0.8 


e2 


110.506127 


12 


e3 


-193.032628 


50 


e 4 


186.624090 


100 




-80.141986 




yt 


1.587 


fixed 


Vu 


-0.82 


fixed 


K c 


0.2216546 


fixed 
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in Table [J to fits including up to the first order. Reasonable fits, as following from the \ 2 
criterion, could only be be obtained by including three more coefficients Cj or Wj. Moreover, 
these coefficients tended to assume much larger values. For this reason, we prefer the fit up 
to second order in u, although the fit up to first order also yields a satisfactory numerical 
representation of the critical energy and specific heat. 

Only the parameters Ct_, O+, Co and c\ y describing the leading few orders of E and C, were 
about the same for both types of fits. It is thus clear that not too much physical significance 
should be given to the subleading and higher-order parameters given in Table (TJ except that 
they provide a numerical description of E and C in the critical region. 

The relative errors in the amplitudes p_ and p + of the second order term in u, as given in 
Table [U are appreciable, and far exceed those of the first-order amplitudes 6_ and b + . For 
this reason we believe that it is not necessary to include a third-order correction, or other 
terms with approximately the same exponent. 



B. Comparison with existing results 



1. Series expansions 

Numerical evaluation of Eq. (I16p allows comparison with results from series expansions. 
The low-temperature series for the energy is provided by Bhanot et al. [15] up to order 25 
in e~ 2K . The specific heat, as obtained by differentiation of this series, is in good agreement 
with Eq. (|T6l) in the interval 0.39 < K < 0.60. The differences, which are shown in Fig. [I], do 
not exceed 10~ 4 . For K > 0.60, outside the range of the least-squares fit, our representation 
of the specific heat with Eq. (|T6|) is no longer accurate and the differences increase sharply. 
The increasing differences for K < 0.30 are due to the truncation of the low-temperature 
series to 25 terms. 

For temperatures above the critical point, a comparison can be made based on the series 



expansion up to order 46 of the free energy as provided by Arisue and Fujiwara [16], with the 
help of Eq. ([3]). The differences with Eq. (fTB"]) are less than 10~ 4 interval 0.15 < K < 0.19, 
as plotted in Fig. [2l For K < 0.15, outside the range of the fit, Eq. ( fTBl) rapidly loses its 
accuracy. The increasing differences for K > 0.19 are due to the truncation of the series. 
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0.008 



0.006 
0.004 
0.002 


0.35 0.4 0.45 0.5 0.55 0.6 0.65 
K 

FIG. 1: Difference Ai = (Clte — Cflt)/^ between the specific heat of the Ising model as obtained 
from the low-temperature series of the energy and from the present least-squares analysis according 
to Eq. (]16p . The difference Ai is at most 10 -4 in the interval 0.39 < K < 0.60. 

0.008 
0.006 

< 

0.004 
0.002 


0.14 0.16 0.18 0.2 
K 

FIG. 2: Difference A^ = (Chte — Cfit)A between the specific heat of the Ising model as obtained 
from the high-temperature series of the energy and from the present least-squares analysis according 
to Eq. (|16p . The difference A^ is at most 10 -4 in the interval 0.15 < K < 0.19. 

2. Amplitude ratios and analytic background 

The fit up to first order in u yielded a universal amplitude ratio a_/a + = 0.540 (5), 
which is to be compared to the result of the fit including the second order of u, which is 
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a_/a+ = 0.532 (7) as follows from the parameter values in Tabled! Based on the consistency 
between these two results, we believe that the latter result a /a+ = 0.532 (7) is reliable. 
This result is close to an estimate 0.541 (14) by Bagnuls et al. 24( from field theory, and to 
the result 0.523 (9) obtained by Liu and Fisher 25] based on series expansions, and slig htly 
smaller than 0.560 (10) as determined from Monte Carlo data by Hasenbusch and Pinn 26] . 

Another universal ratio that can be constructed from the results in Table [J concerns the 
corrections-to-scaling amplitudes. The data in the table suggest 6-/6+ = —0.34 (9), which 
differs considerably from 6-/6+ = —0.96 (25) as obtained by Bagnuls et al. [24j (note the 
sign difference with respect to the notation used by Bagnuls et al., which is related to the 
factor d\t\/dK in our Eq. (1161) ). The sign of this amplitude ratio is in agreement with the 



conclusions of Liu and Fisher 



27J. 



As noted in Sec. IIV A\ there may be corrections to scaling governed by an irrelevant field 
u' with exponent y u i ~ 2y u , and thus indistinguishable from contributions in second order of 
u. It is thus possible that the amplitudes p + and p_ as given in Tabled contain contributions 
due to the field u' . Therefore, the resulting ratio P-/p+ = 0.61 (24) may not qualify as a 
universal amplitude ratio. 

Our result for the critical energy, = 0.990604 (4), can be compared with results 
obtained from series analysis. It is slightly smaller than the result eo = 0.99218 (15) obtained 



by Sykes et al. 28], slightly larger than eo = 0.9902 (1) found by Liu and Fisher [251 ]. and in 
agreement with e = 0.991 (1) found by Butera and Comi [29|. Our result is also consistent 
with the Monte Carlo estimates eo = 0.990 (4) due to Jensen and Mouritsen [30(, and 
e = 0.9904 (8) due to Hasenbusch and Pinn 



3. Comparison with experimental results for Rb^CoCk 



As implied in the Introduction, the magnetic Co 2+ ions in rubidium cobalt chloride 
assume a spin-1/2 Ising character. This has been experimentally confirmed 3jj in the 
related compound CS3C0CI5. The magnetic moments are aligned along the c direction of the 
tetragonal crystal structure. The Co 2+ ions are arranged in a simple Bravais lattice, with 
equivalent positions j^]. Furthermore, electron-spin resonance results 33] for CS3C0CI5 
showed that the exchange interaction with the two nearest neighbors in the crystallographic 
c direction has the same magnitude as that with the four nearest neighbors in the aa plane, 
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so that one may expect that the theoretical results for the simple-cubic Ising model are 
applicable. Specific-heat and magnetic susceptibility measurements js] on RbsCoCls showed 
that a phase transition to an antiferromagnetic phase occurs at T c = 1.14 K. It was indeed 
found that the specific heat (which does not depend on the sign of K) did agree with 
the theoretical predictions available at that time. These predictions were based on series 

35I ]. and on the assumption that the specific- heat 



expansions due to Baker [34| and Sykes 
exponent a = 0. In view of later results for the specific- heat exponent, as well as the effect 
of Wegner's correction the comparison made in Ref. ^ may thus not be considered 



as entirely satisfactory. In Fig. [3] we show the experimental data together with Eq. (TTBT) . 
as well as results from the low- and high-temperature series. This comparison with the 
experimental data, which involves only one adjustable parameter, the critical temperature, 
shows that the specific heat of RbsCoCls agrees reasonably well with the predictions for the 
simple-cubic Ising model. The data in Fig. [3] suggest small deviations at low as well as at high 
temperatures, but there the specific heat becomes very small, so that the experimental error 
margins, which include the uncertainty of heat capacity of the empty apparatus, become 
appreciable. A comparison of the experimental data listed in Ref. |9| with the results from 
Eq. ffTBT) show that deviations up to a few percent occur also in the range 0.9 < K/K c < 1.2. 
But these deviations do not display an obvious systematical trend, and may possibly be 
attributed to the fact that the measured heat capacity is, near criticality, the result of 
integration of a highly nonlinear function over a nonzero temperature range. 

It thus seems that new experiments on RbsCoCls are needed to firmly establish deviations 
with respect to the predictions for the simple-cubic Ising model. Such deviations would be 
a logical consequence of the tetragonal symmetry of RbsCoCls, which implies that there 
is no reason why the coupling in the c-direction should be precisely equal to that in the a 
direction. Also the presence of interactions with further neighbor spins, which include small 
magnetic dipole-dipole interactions, should lead to deviations. 



C. Conclusion 



The formula Eq. (TIB]) , supplemented by Eqs. ( fT7l) . ( fl8l) . and (TP2l) and by the parameter 
values in Table (U describes the specific heat of the three-dimensional Ising model in the 
interval 0.15 < J/kT < 0.60. Comparisons with low- and high-temperature series expansions 
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FIG. 3: Specific heat of the Ising model on the simple-cubic lattice. Logarithmic scales are used 
because of the large variation of the specific heat with temperature. The data points are existing 
experimental results [9|] for RbsCoCls. The full line represents the scaling form Eq. (|16|) with the 
parameters defined in Table HI The dashed lines at the lower left and right are obtained from low- 
and high-temperature series expansions [15I . of the the free energy. 



yield satisfactory agreement in the intervals 0.15 < J/kT < 0.19 and 0.39 < J/kT < 0.60 
respectively. The differences between Eq. (fl~6l) and the results from series expansions are 
at most 10~ 4 in the mentioned intervals. These differences are smaller than the statistical 
errors in the Monte Carlo results on which Table [J is based, as may be expected since the 
number of 100 data points far exceeds the number of 16 free parameters in the fit formula, 
so that in effect averaging occurs. Since Eq. ([161) continues to satisfactorily describe the 
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extrapolated Monte Carlo data until a distance \K — K c \ m 0.005 from the critical point, we 
conclude that the error margin in Eq. (fT6j) does not exceed that of the Monte Carlo data, 
i.e., it will be limited to at most a few times 10~ 3 at \K — K c \ > 0.005. Larger uncertainties 
are expected for \K — K c \ < 0.005 because of the error margins in the critical amplitudes, 
exponents and temperature. Taking into account these numerical uncertainties, Eq. HTB]) 
can be used in the interval 0.15 < J/kT < 0.60 for comparison with experiments on systems 
that are described by the simple-cubic Ising Hamiltonian. 

In addition, our results show that Monte Carlo simulations can be used to determine the 
universal leading amplitude ratios a_/a + and even the nonasymptotic ratio Thus 
far, the correction amplitudes have been studied by means of series analysis, field theory, 



and crossover scaling, see e.g., Refs. 
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TABLE II: Extrapolated values of the dimensionless energy density of the infinite, simple-cubic 
Ising model as a function of the coupling K. Estimated error bounds are included. 



K 


-E{K) 


error 


K 


-E{K) 


error 


0.12 


0.382236 


0.000020 


0.13 


0.419100 


0.000109 


0.14 


0.457696 


0.000007 


0.15 


0.498271 


0.000098 


0.16 


0.541261 


0.000008 


0.166 


0.568535 


0.000448 


0.17 


0.587433 


0.000007 


0.172 


0.597108 


0.000047 


0.178 


0.627251 


0.000049 


0.18 


0.637719 


0.000009 


0.184 


0.659300 


0.000052 


0.19 


0.693640 


0.000008 


0.195 


0.724504 


0.000007 


0.196 


0.730987 


0.000060 


0.2 


0.757945 


0.000005 


0.202 


0.77224 


0.000011 


0.205 


0.794807 


0.000016 


0.21 


0.836533 


0.000053 


0.213 


0.865007 


0.002836 


0.22165460 


0.990604 


0.000004 


0.224 


1.135886 


0.000083 


0.225 


1.184077 


0.000113 


0.226 


1.228851 


0.002835 


0.227 


1.270902 


0.001567 


0.228 


1.310761 


0.000069 


0.229 


1.348806 


0.000084 


0.23 


1.385036 


0.000010 


0.232 


1.453408 


0.000128 


0.235 


1.547095 


0.000082 


0.24 


1.684411 


0.000006 


0.25 


1.908377 


0.000006 


0.26 


2.084422 


0.000011 


0.268 


2.200165 


0.000117 


0.27 


2.226207 


0.000025 


0.28 


2.342331 


0.000013 


0.29 


2.438392 


0.000020 


0.3 


2.518570 


0.000013 


0.31 


2.585908 


0.000017 


0.32 


2.643039 


0.000137 


0.33 


2.691193 


0.000121 


0.34 


2.732442 


0.000779 


0.35 


2.767639 


0.000008 


0.36 


2.797886 


0.000224 


0.37 


2.823926 


0.000119 


0.38 


2.846340 


0.000055 


0.39 


2.865798 


0.000070 


0.4 


2.882622 


0.000005 


0.42 


2.909917 


0.000030 


0.44 


2.930623 


0.000034 


0.45 


2.939050 


0.000003 


0.46 


2.946380 


0.000022 


0.48 


2.958427 


0.000035 


0.5 


2.967777 


0.000002 


0.6 


2.990703 


0.000001 


0.65 


2.994958 


0.000001 


0.7 


2.997255 


0.000001 
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TABLE III: Extrapolated values of the dimensionless specific heat of the infinite, simple-cubic Ising 
model as a function of the coupling K. Estimated error bounds are included. 



K 


C(K)/k 


error 


K 


C(K)/k 


error 


0.12 


0.05212 


0.00004 


0.14 





07736 


0.00005 


0.15 


0.09382 


0.00013 


0.16 





11381 


0.00004 


0.166 


0.12812 


0.00008 


0.17 





13882 


0.00005 


0.172 


0.14458 


0.00010 


0.178 





16391 


0.00012 


0.18 


0.17094 


0.00012 


0.19 





21445 


0.00035 


0.195 


0.24361 


0.00031 


0.2 





27950 


0.00110 


0.205 


0.32850 


0.00093 


0.208 





36500 


0.00100 


0.21 


0.39573 


0.00114 


0.214 





48216 


0.00165 


0.215 


0.51335 


0.00261 


0.216 





54767 


0.00157 


0.217 


0.58500 


0.00200 


0.227 


2 


11009 


0.00675 


0.228 


2.01738 


0.00323 


0.229 


1 


95000 


0.01000 


0.23 


1.87829 


0.00300 


0.232 


1 


76800 


0.00150 


0.235 


1.64500 


0.01000 


0.238 


1 


52800 


0.00300 


0.24 


1.46600 


0.00200 


0.25 


1 


23138 


0.00041 


0.26 


1.06294 


0.00031 


0.27 





93268 


0.00200 


0.28 


0.82520 


0.00060 


0.29 





73597 


0.00094 


0.3 


0.66021 


0.00072 


0.31 





59383 


0.00085 


0.35 


0.39923 


0.00039 


0.36 





36325 


0.00036 


0.37 


0.33045 


0.00024 


0.38 





30110 


0.00020 


0.39 


0.27460 


0.00020 


0.4 





25058 


0.00025 


0.42 


0.20890 


0.00020 


0.44 





17428 


0.00013 


0.45 


0.15910 


0.00020 


0.46 





14567 


0.00012 


0.48 


0.12158 


0.00010 


0.5 





10160 


0.00010 


0.55 


0.06470 


0.00015 


0.6 





04109 


0.00010 


0.65 


0.02593 


0.00008 


0.7 





01632 


0.00004 
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